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ABSTRACT 

A series of stellar models of spectral type G is computed to study the rotation laws resulting from mean-field equations. The rotation 
laws of the slowly rotating Sun, the fast rotating MOST stars e Eri and Cet and the rapid rotators R58 and LQ Lup can easily be 
reproduced. We also find that differences in the depth of the convection zone cause large differences in the surface rotation law and 
that the extreme surface shear of HD 171488 can only be explained with a artificially shallow convection layer 
We also check the thermal wind equilibrium in fast-rotating G dwarfs and find that the polar subrotation (d£l/dz < 0) is due to the 
barocline efi"ect and that the equatorial superrotation (dfi/dr > 0) is due to the A effect as part of the Reynolds stresses. In the bulk 
of the convection zones where the meridional flow is slow and smooth the thermal wind equilibrium actually holds between the 
centrifugal and the pressure forces. It does not hold, however, in the bounding shear layers including the equatorial region where the 
Reynolds stresses dominate. 
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1. Introduction 

Many stars show signs of differential rotation. The solar equa- 
tor rotates with a shorter period than the polar caps. The dif- 
ference of 132 nHz between the rotation periods found by 
[Ulrich et al. (1988)1 corresponds to a difference of 0.07 rad/day 
between the respective angular velocities or a lapping time of 
88 days. Helioseismology has found that this pattern persists 
throughout the whole convection zone but not in the radia- 
tive zone below ( [Thompson et al. 2003] l. Stellar differential ro- 
tation can be inferred from the light curves of rotating spotted 
stars (see, e.g. Henry et al. 1995; Messina & Guinan 2003), 
from monitoring magnetic activity (Donahue et al. 1996), spec- 
troscopically (Reiners & Sch mitt 2003 ), or by Doppler imaging 
'. jBarnes et al. ioOOI IEonati "et al. 2000l l. 

While several studies have found a systematic dependence of 
, the surface differential rotation on the rotation rate, no such de- 
' pendence is found when the samples are combined (Hall 1991; 
I IBarnes et al. 20051 ). Moreover, measurements of surface differ- 
, ential rotation of the rapidly rotating K dwarfs PZ Tel and AB 
■ Dor with the Doppler imaging technique show that stars rotating 
much faster than the Sun show very similar surface shear values 
- as first predicted by Kitchatinov & Riidiger (1999). PZ Tel is 
a young K dwarf with a rotation period of 0.95 days. Its surface 
differential rotation 6S1 = 0.075 rad/day is remarkably close to 
that of the Sun dBarnes et al. 2000). AB Dor is a rapidly rotating 
KO dwarf with a rotation period of 0.5 1 days. Its surface differ- 
ential rotation was found to vary with time between 0.09 and 
0.05 rad/day (ICoUier Cameron & Donati 20021 ). 

Barnes et al. (2005) proposed a dependence on the effective 
temperature (albeit with large scatter) and found the power law 



where 6Q. is the difference between the angular velocities at the 
equator and the polar capsQ. The power law ([TJ was confirmed 
by the findings using spectroscopic methods (Reiners 2006). The 
light curves of the stars e Eri and /<•' Cet recorded by the MOST 
satellite dCroll et al. 20061 IWalker et al. 200"7] ) both indicate with 
- 0.062 and SO. ^ 0.064 a very similar equator-pole differ- 
ence of the rotation rate as the Sun. Also the value of 0. 1 1 for 
the young G star CoRoT-2a (Frohlic h etal. 20091 ) well fits the 
common picture of a rotation-independent surface shear for G 
stars. 

It is challenged, however, by recent observations of the 
young G dwarfs LQ Lup, R58 and HD 17 1488. Marsden et 
al. (2005) report a surface differential rotation of 0.025 + 0.015 
rad/day for R58 while |Jeffers & Donati (2009b)l find the much 
larger value of 0.138 + 0.011 rad/day. Donati et al. (2000) find 
a surface differential rotation 6Q = 0.12 + 0.02 rad/day for LQ 
Lup. Jeffers & Donati (20 09a) determined the surface rotation of 
HD 171488 using the Zeeman-Doppler imaging technique and 
found a very strong surface differential rotation of 0.5 rad/day. 
Marsden et al. (2006) report a smaller but still large value of 
0.402 + 0.044 rad/day. Huber et al. (2009) found no evidence for 
differential rotation at all but could not rule it out either Jeffers 
et al. (2010) confirmed the findings of Marsden et al. (2006). 

The values found for LQ Lup and R58 are in line with the 
Barnes et al. picture but the large values found for HD 171488 
are not. The three G dwarfs are similar in their ages, effective 
temperatures and radii yet HD 171488 shows a much stronger 
differential rotation than the other two stars. The studies men- 
tioned have focused on rotation rate and effective temperature as 
the properties determining the surface differential rotation. As 
the stars are very similar by their stellar parameters such as age 

' the strength of the differential rotation is also expressed in terms of 
the lapping time between the equator and the poles, Pi^p = 2n/SQ. 
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and effective temperature, we ask if there could be a difference 
in their internal structure that would cause the observed differ- 
ence in their surface rotation. All three stars have just reached 
the zero age main sequence or are approaching it. Given the rapid 
retreat of the convection zone in the final part of the pre-main se- 
quence evolution and the uncertainty of the stellar age we ask if 
the strong differential rotation of HD 1714888 can be explained 
by a different depth of its outer convection zone. 

In the following we compute model convection zones of dif- 
ferent depth and their large-scale gas motions, i.e. rotation and 
meridional flow. The models are based on the mean field formu- 
lation of fluid dynamics which has been very successful for the 
Sun, where the models reproduce the surface rotation, the inter- 
nal in the convection zone, and the surface meridional flow very 
well (Kitchatinov & Rudiger (1999, KR99), Kuker & Stix (2001, 
KSOl)). Models have been constructed for a variety of stars with 
spectral types from M to F (Kiiker & Rudiger 2008). 

A new scheme allows the computation of stellar rotation 
laws and meridional flow patterns based on a mean-field model 
of the large-scale flows in stellar convection zones also for fast 
rotation rates when narrow boundary layers exist. It assumes 
strict spherical symmetry for the basic stratification, ignoring 
any flattening that might occur for very fast rotation. However, 
flie impact of rotation on the thermal structure can be taken into 
account by including a gravity darkening term in the heat trans- 
port equation so that the model remains applicable even for mod- 
erately flattened stars. 

2. Theory of stellar rotation laws 

While thermal convection is driven by the stratification of the 
star, the convective gas motions carry angular momentum as well 
as heat. Momentum transport by turbulent motions is known as 
Reynolds stress and can be described as a turbulent viscosity in 
case of the simplest shear flows. In a rotating, stratified convec- 
tion zone the Reynolds stress is anisotropic and its azimuthal 
components have a contribution proportional to the angular ve- 
locity, Q, itself rather than its gradient. This form of Reynolds 
stress (with 'A effect') is not compatible with solid body rota- 
tion. 



2.1. Basic equations 

Our model consists of a ID background model which assumes 
hydrostatic equilibrium and spherical symmetry and a system 
of partial differential equations that describe the convective heat 
flux, the transport of angular momentum, and the meridional 
flow, assuming axisymmetry. The equation for the convective 
heat transport then reads 

V . (F™"^ + F"'!) - pTu ■ V5 = 0, (2) 

where /i'<=°'™ and f are the convective and radiative heat flux, 
respectively, p the mass density, u the mean gas velocity, T the 
gas temperature. S is the specific entropy, 

5=C.,log(^j, (3) 

where p is the gas pressure and C,. the specific heat at constant 
volume and y the adiabate index. In spherical polar coordinates 
the Reynolds equation can be rewritten as a system of two partial 
differential equations. The azimuthal component of the Reynolds 
equation expresses the conservation of angular momentum. 



where t is the angular momentum flux vector 
tj - r sin 6(pr sin OQ.u'^ + pQi^), 



(5) 



with the two transporters of angular momentum i) the meridional 
flow (m™) and ii) the zonal flux of the turbulent angular momen- 
tum (Qi^). 

The equation for the meridional flow is derived from the 
Reynolds equation by taking the azimuthal component of its 
curl: 



P 



do? 1 

-1- r sin 0—— + -r(Vp X Vp)^ -h . . . = 
oz 



(6) 



where Rij — -pQij is the Reynolds stress and d/dz = cos 6d/dr- 
sin 68/ do is the derivative along the axis of rotation. For very fast 
rotation the second term on the LHS of Eq. ^ dominates and the 
rotation rate is constant along cylindrical surfaces, as stated by 
the Taylor-Proudman theorem. 

2.2. Transport coefficients 

If spherical polar coordinates are chosen, the A effect appears in 
two components of the Reynolds stress: 



Qr<^ = e7 + yvtsin0Q 
Qe^ = e™=+//y,cos0Q 



(7) 
(8) 



Here, 2^^^'^ and Qg'^'^ contain only first order derivatives of Q. 
with respect to r and 6 and therefore vanish for uniform rota- 
tion. The coefficients V and H refer to the vertical and horizontal 
part of the A effect while Vt represents the eddy viscosity in the 
convection zone. As the corresponding terms contain the angu- 
lar velocity Q. itself, they do not vanish for rigid rotation. Thus, 
rigid rotation is not stress-free. 

For slow rotation, the convective heat flux is proportional to 
the gradient of the specific entropy: 



-pTXt 



dS_ 

dxi' 



(9) 



where ;^'t the turbulent heat conductivity. The turbulent heat con- 
ductivity is determined by the convection velocity and the 
convective turnover time Tr'. 



1 2 

Using standard mixing-length theory we find 

2_ ligds 



4Cp dr 



and 



(10) 



(11) 



(12) 



where is the mixing length. Equation ^ describes a strictly 
radial, diffusive heat flux. In a rotating convection zone the heat 
flux is not strictly aligned with the entropy gradient, i.e. 



(13) 



V ■ ^ = 0, 



(4) 



(Rudiger 1989, [Kitchatinov et al. 199"4l ), where the dimension- 
less coefficients 0,j are functions of the Coriolis number Q* — 
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2Tcfi fulfilling lO/jl < 1. The flux vector is then tilted towards 
the rotation axis. 

The radiative heat flux is given by 



7rad 



l6o-T^ViT 
3kp 



(14) 



where cr is the Stefan-Boltzmann constant and k the opacity. 



2.3. Background model 

We assume that the star is essentially in hydrostatic equilibrium 
and that all gas motions constitute small perturbations of that 
equilibrium which do not change the structure of the star. We 
particularly assume that no net gas transport occurs: 



V ■ (pS) = 



(15) 



Our model convection zone is a spherical layer with adiabatic 
stratification that is heated from below and cooled at the surface. 
We assume an ideal gas with the equations of state. 



P = -pT = Cp pT 

fi y 



(16) 



(p is the gas pressure, fi the gas constant, Cp the heat capacity 
for constant pressure) and hydrostatic equilibrium. 



dp 



-gP- 



From the law of gravity, 
GM{r) 



(17) 



(18) 



(G gravity constant, M(r) the mass contained in the sphere of 
radius r), we derive 



dr r 



(19) 



For adiabatic stratification, pp ^ - const., the density can be 
expressed in terms of the temperature. 



P=Pe - 



i/(y-i) 



(20) 



where pe and are the values of density and temperature at a 
reference radius r^. The condition of hydrostatic equilibrium can 
then be rewritten in terms of the temperature, i.e. 



(21) 



Equations (fT9l l. (l20l i. and (ISTT i form a system that can be inte- 
grated from the reference radius, where the values of density, 
gravity, and temperature are known from the full stellar evolu- 
tion model. 



2.4. Boundary conditions and numerical method 

As boundary conditions we assume stress-free and impenetrable 
boundaries for the gas motion and an imposed radial heat flux. 
At the lower boundary the heat flux is constant; corresponding 
to the stellar luminosity for the heat transport. 



(22) 



where L is the stellar luminosity and r\, the radius at which the 
boundary is located. As the radiative heat flux is prescribed and 
we define the bottom of the convection zone as the point where 
the radiative heat flux equals the total heat flux this condition is 
equivalent to imposing zero convective heat flux. 

At the upper boundary, we assume that the convective 
heat flux equals the radiative flux is converted into radiation 
dGilman & Glat z maier 1981,) , which is expressed through the 
linearized condition 



L S 



(23) 



where rx is the radius of that boundary. The boundary conditions 
on the rotation axis are implied by the requirements for axisym- 
metry. 

To solve the above system an expansion in terms of spherical 
harmonics is used for the dependence on latitude. For the specific 
entropy, the angular velocity and the stream function ^ of the 
meridional flow. 



1 5iA 
pr^ sm 6 o8 

we have 

N 

S{r,e) = J]5„(r)P2(«-i)(cos0), 



1 dilr 
pr sin 6 dr ' 



(24) 



n=l 
N 



P2„_,(COS0) 



Q(r,0) = ya,(r)- . , 

n=l 
N 

ilj{r,6) = J])A«('-)^L(cos0)sin0. 



(25) 



In this equation. Pi and P, are the normalized standard and ad- 
joint Legendre polynomials. The expansions are convenient be- 
cause the functions of colatitude on the right hand side are the 
eigenfunctions of the angular parts of the corresponding diffu- 
sion operators. In using only Legendre polynomials of only even 
or odd degree in a certain expansion we assume symmetry with 
respect to the equator. The expansion ( |25] | usually converges fast 
such that = 5 is sufficient for cases like the solar rotation. 
Faster rotating stars require values up to 20. 

Applying the expansions ( l25l l to the system of partial differ- 
ential equations (|2|i, Q, and (|6]l transforms the equations into a 
system of ordinary nonlinear differential equations with the in- 
dependent variable r. We further reduce the equations to the sys- 
tem of first-order differential equations by introducing new de- 
pendent variables. In particular, the Reynolds stress components 
Rro and R^^ are convenient new dependent variables. The result- 
ing system of first-order differential equations is solved by the 
standard relaxation method as described by Press et al. (1992). 

The rapidly rotating stars have very thin boundary layers 
near the top and bottom of their convection zones. To resolve the 
layers, nonuniform numerical grid (zeros of Chebyshev polyno- 
mials) is applied 



1 / ^ , , / '-3/2 
ri = 2 I ''t + '"b - (n - '"b) cos U ^ _2 



(26) 



2 < / < n - 1, 



ri = rb, r„ = n , 



where n is the number of radial grid points. The grid has fine 
spacing near the boundaries. To solve for the solar rotation law, 
a value of n as small as 20 is sufficient, but we usually use higher 
resolutions (e.g. n - 200), especially for fast rotation. 
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3. The Sun 

We first simulate the solar rotation law and compare the re- 
sults to those from earlier models by Kitchatinov & Riidiger 
(1999) and Kuker & Stix (2001) solving the same set 
of equations with different numerical methods. KR99 used 
a simple Kramers' law for the opacity while KSOl used 
the opacity from Ahrens et al. (1992). Now an analyti- 
cal opacity law after Stellingwerf (1975)1 as distributed with 
[Hansen & Kawaler (1994) is used. A more elaborate treatment 
is possible but not likely to improve the model as we do not 
treat the uppermost layers of the convection zone and the atmo- 
sphere. As the solar convection zone is not fully ionized up to the 
photosphere, the adiabatic gradient, V^d = (dlogT /dlogP)^ , 
and hence the specific heat capacity are not constant. While 
Vad - 0.4 holds throughout most of the convection zone, it drops 
to values as low as 0. 1 in the upper 35,000 km. As we work with 
constant values we either have to exclude the uppermost layer or 
adjust Vad to reproduce the depth of the convection zone from 
the stellar evolution model. 




0.19 0.21 0.23 0.24 0.26 FRACTIONAL RADIUS 




20 40 60 80 
LATITUDE 



Fig. 1. The solar differential rotation and meridional flow as 
computed with the new scheme. Top left: Surface rotation rate 
vs. co-latitude. Top right: rotation rate as function of radius at 0, 
15, 30, 45, 60, 75, and 90° latitude, respectively, from top to bot- 
tom. Bottom left: Streamlines of the meridional flow. The flow 
is counter-clockwise, i.e. directed towards the pole at the surface 
and towards the equator at the bottom. Bottom right: Flow speed 
as function of latitude at the top (solid blue) and bottom (dashed 
red) of the convection zone. Positive values indicate a flow away 
from the north pole. 

The top panels in Fig. [T| show the rotation pattern in the so- 
lar convection zone. It agrees well with those computed with the 
KR99 codes and the KSOl scheme with the modifications de- 
scribed in Bonanno et al. (2007). The surface rotation is fastest 
at the equator with a difference 

(5n = 0.07 rad/day (27) 



BOTTOM OF CZ TOP OF CZ 




20 40 60 SO 20 40 50 80 
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Fig. 2. Temperature deviation as function of the latitude at the 
bottom (left) and top (right) of the solar convection zone. The 
polar axis is always warmer than the equator The differences, 
however, are small. 

between the rotation rates of the equator and the polar caps, cor- 
responding to a relative shear, 5Q/Qeq, of 29 percent. The vari- 
ation with radius at midlatitudes is weak. The angular velocity 
decreases with radius at high latitudes while at the equator it 
slightly increases with radius ('superrotation') in the deep con- 
vection zone. In the surface layer there is a negative gradient of 
Q ('subrotation') at all latitudes. Like the rotation profiles com- 
puted with the KR99 and KSOl schemes, this rotation profile is 
in excellent agreement with the findings of helioseismology. 

The bottom panels of Fig.[T]show the meridional flow. There 
is one cell per hemisphere with the surface flow directed towards 
the poles and the return flow at the bottom of the convection 
zone. The amplitude of this 'counter-clockwise' flow is 14 m/s 
at the (model) surface and 6 m/s at the bottom. The difference 
between the flow speeds at the top and bottom respectively is 
smaller than might be expected from the density stratification 
and the requirement of mass conservation but the concentration 
of the return flow to a thin layer allows a relatively fast flow 
despite the larger mass density. 

Figure |2] shows the quantity 

6T^(S-So)^, (28) 

at the top and the bottom of the solar convection zone with S o 
the specific entropy at the bottom of the convection zone at the 
equator. This choice of 5 implies negative values of 6T at the 
top of the convection zone. In general, the polar axis is warmer 
than the more equatorward parts of the fluid. Unfortunately, the 
smallness of the temperature difference does not allow empirical 
confirmations. Though the deviation from the adiabatic stratifi- 
cation is only small, the resulting barocline term has a profound 
impact on the large-scale meridional flow and the differential ro- 
tation (see below). 

3.1. Fast-rotating Sun 

To study the impact of the basic rotation we compute the rota- 
tion law of a hypothetical fast-rotating Sun with the P = 1.33 
day rotation period of HD 171488. Figure [3] shows the resulting 
rotation pattern. The isocontours are cylinder-shaped in accord 
to the Taylor-Proudman theorem. Consequently, the radial pro- 
files in the right diagram show an increase of the rotation rate 
with increasing radius at low latitudes. Unlike the case of the 
real Sun, there is a pronounced increase in the upper part of 
the convection zone at the equator At both radial boundaries 
there are pronounced boundary layers with huge gradients of 
the rotation rate which are caused by the stress-free boundary 
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conditions. The surface rotation is fastest at the equator and de- 
creases monotonously towards the polar caps but the slope of 
the decline has a minimum at mid-latitudes. With a value of 
0.08 rad/day instead of 0.07 rad/day the horizontal shear is only 
slightly stronger than for the Sun despite a factor of 20 between 
the average rotation rates. 

The meridional flow only has one flow cell per hemisphere 
with the surface flow directed towards the poles. The amplitudes 
are 34 m/s at the surface and 17 m/s at the bottom of the con- 
vection zone. As for the real Sun with its much slower rotation, 
the return flow is half as fast as the surface flow but even more 
concentrated at the bottom of the convection zone. The flow is 
fastest at mid-latitudes. 





4.67 4.69 4.71 4.73 4.75 



0.750.800.850.900.95 
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0.69 For this model star we find a surface differential rota- 
tion 6Q. = 0.51 rad/day or A: = 0.10. This is in agreement with 
the observed value. 

For A-' Cet we use solar mass and metalicity and an age of 
1 Gyr The model star has an effective temperature of 5677 K, 
0.915 Rq and 0.78 Lq. The resulting rotation law shows a sur- 
face shear 6Q. - 0.077 rad/day, or ^ = 0.12. This is slightly 
larger than the observed value. An earlier model from the same 
evolutionary track, with a stellar age of 167 Myr, an effective 
temperature of 5649 K, a radius of 0.9 Rq, and an luminosity of 
0.74 Lo yields k = 0.11, or 6Q. = 0.072 rad/day. This is still 
not in perfect agreement with the observed value but reasonably 
close. 

Possible reasons for the remaining discrepancy are the stellar 
model, which might not be a sufficiently accurate representation 
of the real star and an underestimation of the total surface shear 
by the sin^ derived from spot rotation as the observed spots do 
not cover the whole range of latitudes from pole to equator. 

4. Numerical experiments 

We carry out several numerical experiments to compare the roles 
played by the two drivers of differential rotation, i.e. the A effect 
and the barocline flow. 



Fig. 3. Rotation of a hypothetical fast-rotating Sun with P 
1.33 d. 



3.2. MOST stars: e Eri and Cet 



- 4.1. Sun without A effect 

Our model includes two effects that are capable to maintain dif- 
ferential rotation, the A effect and the barocline flow as due to a 
horizontal temperature gradient. This can be seen when the baro- 
cUne term in Eq. (|6]l is rewritten in terms of the temperature. 



e Eri and /c' Cet are active dwarf stars for which surface differen- 
tial rotation has been derived from light curves recorded by the 
MOST satellite Croll et al. (2006) derived a rotation law of the 
form 



eq 



1 -ksin^B 



(29) 



for the K2 V star e Eri, where P(B) is the rotation period at the 
latitude B, and Pgq the rotation period at the equator (B -0). The 
parameter k gives the shear of the rotation law. The general rule 
derived from the theory of the stellar rotation laws of Kitchatinov 
&Rudiger(1999) 



eq 



100 days 



(30) 



is perfectly fulfilled by these stars. 

An equatorial rotation period Peq = 11-2 days and a value 
k - 0.11^0 02 for the differential rotation has been reported for 
6 Eri corresponding with 6Q. = 0.062 rad/day. For the G5 V 
star Cet Walker et al. (2007) found Peq = 8.77 days and k = 
0.09+° °°^ con-esponding with 5Q. = 0.064 rad/day. 

Kitchatinov & Olemskoy (2010) computed rotation laws for 
model stars with 0.8 and 1.0 solar masses for e Eri and /c' Cet. 
The results were k - 0.127 for e Eri and k = 0.13 for a:' Cet. 

Here the calculations for these stars are based on im- 
proved stellar models computed with the MESA/STARS code 
jPaxtonet al. 2010b . A star with = O.85M0, Z = 0.02 and 
an age of 0.44 Gyr serves as model for e Eri. It has an effective 
temperature of 5076 K, a radius of 0.76 Rq, and a luminosity of 
0.34 Lq. The bottom of the outer convection zone is located at 



4(Vp X Vp)^ 
P 



j_d6T 
rT 80 ' 



(31) 



The barocline term can be a powerful driver of meridional flows, 
which in turn can drive differential rotation. In our theory the 
latitudinal temperature profile is due to the anisotropic heat- 
conductivity tensor in the presence of rotation. Always the polar 
axis is slightly warmer than the equator but with an unobservable 
temperature excess. A positive pole-equator temperature differ- 
ence will drive a clockwise flow. Its angular momentum transport 
leads to an accelerated rotation at low latitudes and to slower ro- 
tation at the polar caps. A barocline flow (also 'thermal wind') 
can therefore maintain differential rotation with solar-type sur- 
face rotation. For an illustration we repeat our computation for 
the Sun with the A terms canceled within the Reynolds stress. 

The resulting rotation and flow patterns are shown in Fig. 
|4] The rotation is indeed solar-type but with SQ. = 0.04 weaker 
than with A effect and the isocontours are distinctly disc-shaped 
at the poles even in the midlatitudes. The equator region shows 
basically no structure. 

The typical superrotation of the deep convection zone be- 
neath the equator known as a result of the helioseismology does 
not occur The reason is simple: if only a (fast) meridional cir- 
culation transports the angular momentum with a pattern sym- 
metric with respect to the equator then the angular momen- 
tum becomes uniform in the equatorial part of the convection 
zone (except, of course, the boundary layers). Hence, there is 
r^Q ~ const independent of the flow direction but in contradic- 
tion to the observationQ. 



- without A effect a superrotation beneath the equator can only be 
due to an anticlockwise flow which is sufficiently slow 
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With amplitudes of 4.7 m/s at the top and 2.1 m/s at the 
bottom of the convection zone the barocline-induced clockwise 
meridional flow is substantially weaker than in the complete 
model - and it goes in the 'wrong' direction. 
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Fig. 4. Solar rotation law and meridional flow without A effect. 
Top: The rotation law. Note the steep negative gradient of the 
angular velocity at the poles and the almost rigid rotation be- 
neath the equator Bottom: The meridional flow is anti-solar: it 
is positive (equatorwards) at the top of the convection zone and 
negative (polwards) at the bottom of the convection zone. 



4.2. Sun without barocline flow 

Next the barocline term is canceled while keeping the A effect 
which maintains the differential rotation together with the merid- 
ional flow caused by the former through the second term on the 
LHS of Eq. (|6]l ('Biermann-Kippenhahn flow', see Kippenhahn 
1963; Kohler 1970). Figure |5] shows the resulting rotation pat- 
tern. The differential rotation is solar-type, i.e. the rotation pe- 
riod is shortest at the equator and longest at the polar caps. With 
6Q. - 0.014 rad/day the surface rotation is much more rigid than 
observed. The isocontour plot shows a distinctly cylinder-shaped 
pattern in the bulk of the convection zone while at the top and 
bottom boundaries the pattern deviates from the cylinder geom- 
etry and the rotation rate falls off with increasing radius at all 
latitudes. The counterclockwise meridional flow has a very sim- 
ilar geometry as in the full model but it is faster with amplitudes 
of 17 and 8 m/s at the top and bottom. 

The subrotation along the polar axis which is typical for the 
solar rotation law and which is produced by baroclinic flows 
(cf. Fig. nil does not exist here. Kohler (1970) and Riidiger et 
al. (1998) have given several numerical examples of this direct 
consequence of the Taylor-Proudman theorem. 




0.24 0.24 0.25 0.25 0.26 FRACTIONAL RADIUS 



Fig. 5. Solar rotation law without baroclinic terms. The structure 
disappears now at the polar axis but it appears at the equatorial 
region. The equatorial plane shows superrotation but the polar 
axis rotates almost rigid. 



5. Young G dwarfs 

5.1. CoRoT-2a 

This G7 V star is a young Sun, i.e. it has solar mass but is much 
younger with an age of 0.5 Gyrs. The star has been observed 
by the CoRoT satellite and found to have a planet with an or- 
bital period of 1.743 days (Alonso et al. 2008). Besides plane- 
tary transits, the light curve shows periodic variation that is most 
easily explained by a rotating, spotted surface with basic rotation 
period of 4.5 days. Spot modeling using circular spots finds an 
excellent fit assuming three spots and a solar-type surface dif- 
ferential rotation with a difference of 0.11 rad/day (Frohlich et 
al. 2009). 




1.32 1.34 1.37 1.39 1.41 FRACTIONAL RADIUS 



Fig. 6. Rotation pattern of CoRoT-2a, a young Sun rotating with 
P - 4.5 days. 



Our model star is based on a model from an evolutionary 
track for the Sun. The age is 0.5 Gyrs, the radius 0.9 solar radii, 
and the luminosity 75 percent of the solar value. The bottom of 
the convection zone is at a fractional radius x - 0.73. 

Figure|6]shows the resulting rotation pattern. The surface dif- 
ferential rotation of 0.09 rad/day reproduces the value observed 
by Frohlich et al. (2009). The rotation pattern is more cylindrical 
than that of the Sun, but less than that of our fast-rotating Sun. 
Similarly, the boundary layers are more pronounced than in the 
Sun but less than for the fast-rotating Sun. The flow pattern is 
similar to that in the solar convection zone. The amplitude is 
slightly larger with a maximum value of 18.6 m/s at the top and 
10.6 m/s at the bottom of the convection zone. 
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5.2. R58, LQ Lup and HD 1 71488 

We now address the recent differential rotation measurements of 
the fast-rotating G dwarfs R58, LQ Lup and HD 171488. R58 
(HD 307938) is a young active G dwarf in the open cluster IC 
2602. It has a photospheric temperature of 5800 K and rotates 
with a period of 0.56 days dMarsden et al. 2005l l. LQ Lup (RX 
J1508. 6-4423) is a post-T Tauri star with an effective temper- 
ature of 5750 + 50 K and a rotation period of 0.31 days, and a 
radius of 1.22 + 0.12 solar radii (Donati et al. 20001 1. correspond- 
ing to a mass of 1.16 + 0.04 solar masses, and an age of 25 + 10 
Myr. The equator-pole Q-difference is 0. 12 rad/day. HD 171488 
(V889 Her) is a young active G dwarf. Strassmeier et al. (2003) 
found a photospheric temperature of 5830 K, a rotation period 
of 1.337 days, a mass of 1.06 + 0.05 solar masses, a radius of 
1.09 + 0.05 solar radii, and an age of 30-50 Myr Mai-sden et 
al. (2006) find a rotation period of 1.313 days, a photospheric 
temperature of 5800K, and a radius of 1.15 + 0.08 solar radii. 

With 0.4-0.5 rad/day the equator-pole O-difference of 
HD171488 is exceptionally large. The A;-value of 0. 10 misses the 
k - 0.013 predicted by Eq. ( l30b by a factor of 7.5. As the corre- 
sponding factors for R58 and LQ Lup are much smaller (factor 
2-3) the following calculations take the case of HD 171488 as 
the main example. As the stars have similar radii and effective 
temperatures we investigate what impact differences in the inter- 
nal structure have on the surface rotation patterns. We study two 
models that mainly differ in the metalicity and, as a consequence, 
the depth of the convection zone. The models were computed 
with the MESA/STARS code. We use the metalicity as control 
parameter for the depth of the convection zone and vary mass 
and mixing length parameter to adjust radius and effective tem- 
perature to values close to those observed for the three young G 
dwarfs. The rotation period is the 1 .33 day period of HD 17 1488. 
Both model stars are 30 Myr old. 

5.2.1. Deep convection zone 

The first model star has the high metalicity of Z = 0.03. A mass 
of 1.11 solar masses and a value of 1 .6 for the mixing length 
parameter lead to an effective temperature of 5685 K and a ra- 
dius of 1.14 solar radii. The bottom of the convection zone is 
located at x - 0.765 and has a temperature of 1.46 x 10^ K. 
The rotation law is shown in the top part of Fig. |2l The iso- 
contour plot shows cylinder-shaped contours similar to those for 
the fast-rotating Sun. The surface rotation law is solar-type with 
a fast-rotating equator. The surface shear is much stronger. We 
find a value of 0.1 1 rad/day, about 1.5 times the solar value. The 
radial profiles show superrotation beneath the equator and sub- 
rotation along the rotation axis. The boundary layers are much 
less pronounced and the surface layer resembles the real Sun, 
i.e. the rotation rate decreases with increasing radius at all lat- 
itudes. The meridional flow is of solar-type (counterclockwise) 
with the surface flow towards the poles and the return flow lo- 
cated at the bottom of the convection zone. The amphtudes are 
23 and 14 m/s, respectively. 




4.64 4.67 4.69 4.72 4.75 FRACTIONAL RADIUS 




4.30 4.43 4.56 4.68 4.81 FRACTIONAL RADIUS 



Fig. 7. Internal rotation of young G dwarfs with rotation period 
of 1.33 days. Top: deep convection zone, 6Q. = 0.11 rad/day at 
the surface . Bottom: shallow convection zone. 6D. = 0.5 rad/day 
at the surface. 

shear of (5Q = 0.5 rad/day. The contour plot shows a disc-shaped 
pattern with flat profiles at low latitudes and a decrease of the 
rotation rate with increasing radius at the polar caps. The radial 
profiles show pronounced boundary layers at the radial bound- 
aries. The flow pattern is mostly solar-type. There is one large 
flow cell with poleward flow at the top and the return flow at the 
bottom of the convection zone. In addition, there is a small cell 
of clockwise flow at low latitudes. The flow speeds at the top 
reach 18 m/s and 2.2 m/s, respectively. The return flow at the 
bottom of the convection zone reaches 3.5 m/s. 

5.2.3. Truncated convection zone 

The models above have (roughly) the same effective temperature 
but differ in mass and structure. To further isolate the effect of the 
depth of the convection zone, we compute a model of a young 
solar mass star and artificially reduce the depth of the convection 
zone by imposing the lower boundary at a fractional radius x = 
0.88 instead of x = 0.78 and compare the rotation pattern with 
that of the full model. In both cases the top boundary is at x = 
0.98. The truncated convection zone thus has half the depth of 
the full model. In both cases the rotation law is solar-type but 
while the full convection zone produces a surface shear of 0.12 
rad/day, the truncated convection zone has 0.32 rad/day. The full 
model has a temperature difference of 24 K at the top and 124 K 
at the bottom. In case of the truncated model the polar cap is 47 
K hotter than the equator at the top and 254 K at the bottom. 



5.2.2. Sliallow convection zone 

Our second model star has a metalicity of 0.01. A mass of 1.08 
solar masses and a mixing length parameter of 1 .0 lead to a ra- 
dius of 1.14 solar radii and an effective temperature of 5750 K. 
The bottom of the convection zone lies at a fractional radius of 
0.89 and a temperature of 6 x 10^ K. The bottom part of Fig.|2] 
shows the resulting rotation law. The surface rotation has a total 



5.2.4. Very fast rotation 

To illustrate the impact of the rotation rate, we repeat the com- 
putation for the shallow convection zone model with rotation pe- 
riod of 0.33 days. Figure [8] shows the resulting rotation pattern. 
The left plot shows solar-type rotation both in the surface differ- 
ential rotation and in the predominantly radial isocontours. This 
is still not cylinder-shaped but closer to the Taylor-Proudman 
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Fig. 8. Differential rotation of the G dwarf with the shallow con- 
vection zone and a rotation period of 0.33 days. Left: isocontours 
of the angular velocity. Right: the rotation profiles at various lat- 
itudes. 

than the disc-shaped pattern we found for a rotation period of 
1.33 days. The surface shear is 0.38 rad/day, the temperature dif- 
ference between polar caps and equator is 1785 K at the bottom 
and 184 K at the top. The maximum flow speeds are 14 m/s at 
the top and 8 m/s at the bottom of the convection zone. 

6. Discussion 

6.1. Shallowvs. deep convection zone 

Our "shallow" and "deep" convection zone models have similar 
radii and effective temperatures but the stars differ in mass and 
metalicity. We have also chosen different values for the mixing- 
length parameter As a result, the depth of the convection zone 
differs between the models. We find very different results for the 
differential rotation of these convection zones. The model with 
a shallow convection zone shows much stronger differential ro- 
tation than the ones with deeper convection zones. Our "shallow 
convection zone" model does indeed reproduce the very large 
value of 0.5 rad/day found by Jeffers & Donati (2009) while the 
"deep convection zone" model with its smaller value of 0.11 is 
in rough agreement with the values observed for LQ Lup and 
R58. 

The result from the shallow convection zone model is in 
line with previous findings on F stars. Extreme values of sur- 
face shear have been observed for F stars by Reiners (2006) and 
found in theoretical models for stars with very shallow convec- 
tion zones by Kiiker & Riidiger (2007). In the latter study, the 
largest values of surface shear are found for the most massive 
stars studied. These stars are not only the hottest, their convec- 
tion zones are also very shallow. 

To achieve the thin convection zone we had to lower the 
mixing-length parameter to a value of 1 .0. This is not only dif- 
ferent from the value used for the "deep" convection zone model 
but much lower than what is usually used in stellar evolution 
models. Hence, the model is probably not realistic and does not 
directly apply to LQ Lup. Not being specialists in the field, we 
leave the question whether a shallow convection zone is possible 
open but we wonder if strong magnetic fields could inhibit the 
convective heat transport in a way that might be mimicked by 
this choice of parameter 

The "truncated convection zone" model avoids the uncertain- 
ties about the "shallow convection zone" model and shows even 
more clearly the impact the depth of the convection zone has on 
the surface rotation. Both the vertical and horizontal temperature 
gradients are steepest for the model with the shallow convection 
zone and flattest for the model with the deep convection zone. At 



the bottom of the convective zone the differences are 750 K and 
100 K for a convection zone of small or large depth, respectively. 
The corresponding values at the top are 55 K and 19 K. 

6.2. Thermal wind equilibrium 

We always find a higher temperature at high latitudes than at 
the equator This is a consequence of the tilt of the convective 
heat transport vector towards the axis of rotation caused by the 
Coriolis force. In spherical polar coordinates the components of 
the heat flux read 



dS dS 

or o6 

dS dS 

or 09 



(32) 
(33) 



The first term in the horizontal component precludes a purely 
radial stratification. Any variation of the specific entropy with 
the radius will cause a horizontal heat flux and thus build up a 
horizontal gradient. This case is profoundly different from that 
of a latitude-dependent but still purely radial heat flux. The latter 
would have a much smaller impact and would result in a much 
weaker differential rotation (Riidi ger et al. 2003| . 

The impact on the maintenance of differential rotation can 
be seen from the equation for the meridional flow. For fast ro- 
tation (i.e. large Taylor number) Reynolds stress and nonlinear 
terms can be neglected and Eq. ^ is reduced to the thermal wind 
equation 



do. 

2rsin0Q — 

dz 



g_d6T_ 
rT do 



0. 



(34) 



It follows that a gradient of the angular velocity along the axis 
of rotation is needed to balance the horizontal temperature gra- 
dient. Hence, a disc-shaped rotation pattern results in the polar 
area (only) due to the action of the baroclinic flow. For warmer 
poles {6T sinks equatorwards) the z-gradient of Q is negative as 
observed. The meridional flow without the baroclinic component 
only yields dO/dz ^ (Kohler 1970, see also Fig.|5]). Hence, the 
empirical finding of negative dQ/dz at the polar axis by helio- 
seismology proves the existence of baroclinic flows in the solar 
convection zone. 

As we have demonstrated by means of Fig.|4]the superrota- 
tion beneath the solar equator is a direct indication for the action 
of the Reynolds stress in the convection zone. All models with- 
out A-effect but with an clockwise (equatorward at the surface) 
circulation lead to dO/dr ; in the bulk of the convection zone 
beneath close to the equator. Such clockwise flows are able to 
accelerate the equator compared with the poles but they cannot 
produce the superrotation beneath the equator. 

Figure |9] shows results from a detailed computation of the 
terms on the LHS of Eq. ( l34l i and their sum vs. the fractional 
radius for a latitude of 45°. The left panel shows the plot for the 
Sun, the right panel that for the deep convection zone G dwarf 
model. The same quantities were computed for the fast Sun and 
the shallow convection zone G dwarf models. 

For the real Sun the two terms cancel in the bulk of the solar 
convection zone where the meridional flow is slow and smooth 
but not close to the boundaries. For the deep convection zone 
the boundary layers are much thinner but also much more pro- 
nounced. These findings indicate that for fast rotation the merid- 
ional flow is mainly driven by the boundary layers while the bulk 
of the convection zone is in thermal wind equilibrium. In the 
boundary layers the meridional circulation is strongly sheared 
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so that the RHS of Eq. (l34l l (which consists on 3"^ derivatives 
of the meridional velocity components multiplied with the eddy 
viscosity) becomes large. 



Drivers of meridional flow 
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Drivers of meridional flow 
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depth of the convection zone has a big impact on stellar rotation. 
Shallow convection zones produce stronger surface shear. 

The presented mean-field theory for G stars naturally ex- 
plains the rotation laws of the Sun, of the MOST stars e Eri and 
a:' Cet, and of such fast-rotating stars like R58 and LQ Lup. The 
discrepancy between these stars and the much stronger differ- 
ential rotation observed for HD 171488 - if real - is hard to 
explain for stars of similar age and spectral type. If, however, 
HD 171488 had a shallower convection zone for some reason, 
its strong surface differential rotation follows immediately. 

During the pre-MS evolution the convection zone retreats 
from the central region and the originally fully convective star 
forms a radiative zone around the core. This change in the depth 
of the convection zone should be reflected by the surface differ- 
ential rotation. 



Fig. 9. The drivers of angular momentum vs. fractional radius at 
45° latitude for the Sun (left, slow rotation) and the deep convec- 
tion zone G dwarf (right, fast rotation). Dashed line: baroclinic 
term. Dash-dotted line: centrifugal term. Solid line: both. Note 
that the balance ( l34b well holds in the bulk of the convection 
zones of fast rotation stars where the meridional circulation is 
slow and smooth. Equation ( |34] | does not hold in the boundary 
areas where the flow is fast and shearing. 

The numerical experiments show that the baroclinic term is 
an important but not the sole driver of the solar differential ro- 
tation. Cancellation of the A effect leaves a surface shear of 
0.04 rad/day while switching off the baroclinic term reduces it 
to 0.014 rad/day. This is not because the A effect is an inefficient 
transporter of angular momentum but because it is balanced by 
the meridional flow. Canceling the meridional flow altogether in- 
creases the surface shear to 0.18 rad/day. This is a well-known 
effect: meridional flows driven by differential rotation act to re- 
duce it ( Riidiger et al. 1998) ). The baroclinic flow has the oppo- 
site direction and builds up the differential rotation along the po- 
lar axis. The observed superrotation beneath the equator, how- 
ever, cannot be produced by meridional circulation but should 
be a direct consequence of the existence of the A effect. 

7. Conclusions 

Stellar differential rotation is driven by Reynolds stress, by 
the centrifugal-induced ('Biermann-Kippenhahn') flow and (for 
stratified density) by the baroclinic flow. In our model both 
meridional circulations have opposite directions. The barocline 
flow becomes important for faster rotation as the convective 
heat-flux deviates from the radial direction. This deviation can 
be interpreted as a tilt of the heat-flux vector towards the rotation 
axis and causes the poles to be slightly warmer than the lower 
latitudes. The baroclinic flow is responsible for the strong neg- 
ative gradient dQ/dz along the rotation axis, and the Reynolds 
stress produces the typical positive dO/dr along the equatorial 
midplane both known from helioseismology. As a result of both 
impacts the isolines of O known for the Sun are almost radial in 
mid-latitudes. 

Studying real stellar models implies varying a bunch of pa- 
rameters as mass, radius, effective temperature and the depth 
of the convection zone are interdependent. Kitchatinov & 
Olemskoy (2010) studied differential rotation along the lower 
ZAMS for fixed rotation rate and found that the surface shear is 
a function of the effective temperature alone. Using some artifi- 
cial models, we find that for the same effective temperature the 
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